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Abstract 



Using standard numerical Monte Carlo lattice methods, we study non- 
universal properties of the phase transition of three-dimensional c/)^ theory 
of a 2-component real field (j) = ^2) with 0(2) symmetry. Specifically, we 
extract the renormalized values of {(fP')/u and r/v? at the phase transition, 
where the continuum action of the theory is / V(/)p+^r(/)^-|-^(/)^]. These 

values have applications to calculating the phase transition temperature of di- 
lute or weakly- interacting Bose gases (both relativistic and non-relativistic) . 
In passing, we also provide perturbative calculations of various 0(a) lattice- 
spacing errors in three-dimensional O(A^) scalar field theory, where a is the 
lattice spacing. 



I. INTRODUCTION 

A long standing problem is how to compute the first correction ATc, due to interactions, 
to the critical temperature for Bose-Einstein condensation of a very dilute non-relativistic 
homogeneous Bose gas in three dimensions. In Ref . [1] , we review how this problem is related 
to three-dimensional 0(2) 0^ field theory [2], present the results of lattice simulations of that 
theory, and so determine AT^.. The purpose of the present work is to provide details of those 
simulations. The lattice results presented here can also be applied to relativistic Bose gases 
[3,4] and to non-relativistic gases in a trapping potential [5]. 

Three-dimensional 0(2) 0^ theory has the continuum action 



S 



2 ' ^ 



where = (0i,02) is a two-component real field and 0^ = 0^ + For fixed u, we will 
vary r to reach the second-order critical point rc{u) of this model. The shift in the critical 
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temperature of a non-relativistic homogeneous single-species Bose gas is given in terms of 
this theory by [2] 

^ = -^A(^\, (1.2) 
Jo oh n 

where m is the boson mass, n is the number density, and 

A{cP\^[{cP\l-[{cP\]^ (1.3) 

represents the difference between the critical-point value of (0^) for the cases of (i) u non- 
zero and (ii) the ideal gas u—0 (with r — > from above). The result for A(0^)c in the 0(2) 
theory (1.1) can only depend on u and so, by dimensional analysis, it must be proportional 
to u. A primary goal will be to discuss in detail the measurement of the numerical constant 
A(0^)c/m from lattice Monte Carlo simulations of the theory. As reported in Ref. [1], our 
result is 



u 



-0.001198 ±0.000017. (1.4) 



The critical value of r is also useful for various theoretical applications, such as de- 
termining corrections to the value of the chemical potential at the transition [5], or to the 
critical temperature of a trapped gas [5], or to the critical temperature of ultrarelativistic 
(j)^ theory at zero chemical potential [4] . The 0^ interaction is associated with an ultraviolet 
(UV) divergence of the three-dimensional theory and so must be renormalized. If one chooses 
the renormalization scale /2 to be of order u then, by dimensional analysis, the renormalized 
value of rc(/i) must be proportional to u^. The precise scheme used to renormalize r, and the 
precise choice of p,, is a matter of convention. In this paper, we will report a measurement 
of the numerical constant Vc/u^ for r defined by dimensional regularization and modified 
minimal subtraction (MS) at a renormalization scale oi p, — u/3. Our result is 



V? 



0.0019201 ±0.0000021. (1.5) 



One can convert to other choices of by the (exact) identity 

r^(/20 _ r^(AZ.) 2 

V? V? ^9(47r)2 //2- ^ 

In section II, we discuss the lattice action we use, its relationship to continuum fields and 
parameters, how we correct for 0{a) lattice spacing errors, and the algorithms we use for 
simulation. Section III details our procedure for finding the transition, based on the method 
of Binder cumulants. In section IV, we present our initial data, show that we have simulated 
moderately large volumes, and then discuss how to analyze the remaining finite-volume 
corrections by making use of the known critical exponents associated with this universality 
class. The corresponding numerical extrapolations of the finite- volume corrections are given 
in section V. Numerical extrapolation of the continuum limit is presented in section VI. A 
table of all our raw data for various size lattices and values of coupling may be found in 
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Appendix A. The derivations of the 0{a) lattice-spacing corrections used in this paper are 
given in Appendix B. The remaining appendices include various discussions of scaling laws 
used in the text, an analytic calculation of results for small lattice volumes, and a critical 
discussion of one of the early simulations, in the hterature, of the Bose- Einstein condensation 
temperature of dilute non-relativistic gases. 



II. LATTICE ACTION, MEASUREMENT, AND ALGORITHM 

The bare lattice Lagrangian has the form 



(2.1) 



where a is the lattice spacing. In an actual simulation, one invariably chooses lattice units 
where a—1, or equivalently works with rescaled fields and parameters — a^/^$, Tq = a^ro, 
and Uq = auQ. For the sake of presentation, however, we will generally avoid specializing to 
lattice units. 

We work on a simple cubic lattice, and will work in cubic volumes LxLxL with periodic 
boundary conditions, corresponding to {L/af sites. The simplest possible implementation 
of the lattice Laplacian, which we call the "unimproved" choice Vy, would be 

Vu*(x) = [*(^ + «i) - 2*(x) + *(x - ai)] , (2.2) 

i 

where the sum is over unit vectors in the three lattice directions: (1,0,0), (0,1,0), and (0,0,1). 
We use instead a standard improvement, which approaches the continuum Laplacian faster 
for smooth fields: 

Vi'*(x) = a"' + 2«i) + + ai) - |$(x) + |$(x - ai) - j^*(x - 2ai)" . 

i 

(2.3) 

The difference can be seen from the Fourier transforms of the operators and Vj, which 



are 



= 0-^5^(2- 2 cos(aA;j)), (2.4) 

i 

= Z] - ^ cos{aki) + ^ cos{2aki)^ , (2.5) 

and have small k limits 

k^^e + 0{a^k% (2.6) 

^2^A;' + 0(a^A;^). (2.7) 

The unimproved Laplacian has O(a^) error while the improved Laplacian has only 0(0"^) 
errors. 
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A. An unimproved calculation of A(0)^ 

One of our tasks is to calculate the continuum ratio A(0^)/m. In the lattice theory (2.1), 
the free field {uq = 0) result for ($^) is 



2/uo=0 



= 2 



keBZ k"^ 



(2.8) 



where the integral is over the Brillouin zone \ki\ < tt/g. Scaling out the dependence on a, 
we define 



./keBZ fc2 



. - , ■ (2-9) 

Ana 7keBZ fc^ 

For the improved Laplacian, we obtain the value of the constant S by numerical integration: 

E ~ 2.75238391120752. (2.10) 
On the lattice, the most straightforward implementation of the ratio A{(f)'^)/u is then 



Ao($^) _ 1 

Uq Uq 



($2 



2E 



Aira 



(2.11) 



This will approach the desired continuum value as ua ^ 0, but the lattice spacing errors at 
small ua will be 0{ua). 

B. Relationship between lattice and continuum fields and parameters 

One of our goals will be to reduce finite lattice spacing errors, so that we can obtain 
better estimates of the continuum limit with a given coarseness of lattice. To eliminate 
errors at a given order in a, one must not only improve the form of the Laplacian but must 
also perform an appropriate calculation of the relationship between lattice and continuum 
parameters. To this end, we will rewrite our bare lattice action in terms of continuum fields 
and parameters {r,u) as 



, (2.12) 



where the renormalizations Z^, Zr, Sr, and 5u depend on r and u. We explain their derivation 
in Appendix B. For continuum r defined by MS renormalization at a scale p,, we find 



Za, — 



2C2 {uaY ,0, 



5r = a 



5u — a ^ 



[ua] 



2E (ua) 2/6 
3 (47r) 9 \ iia J {Att)^ 



+ 0{{uaf) 



5{ (ua) 



3 (47r) 



32Ci ^ 2^ t""*^' 



9 



(47r) 



0{{ua)^) 



(2.13a) 
(2.13b) 
(2.13c) 
(2.13d) 
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where we have introduced several new numerical constants, given by various integrals in 
lattice perturbation theory, whose values (for the action with the improved Laplacian) are 
approximately 

C ~ -0.083647053040968, (2.14a) 

Ci ~ 0.0550612, (2.14b) 

C2 ~ 0.0334416, (2.14c) 

C3 ~ -0.86147916. (2.14d) 

One needs to similarly match the operator 0^ whose expectation is taken in determining 
A(0^). In Appendix B, we discuss the relationship between the continuum and lattice 
operators and show that the continuum result for A(0^) is 

A(02) = Zr{(t^\,, - + 0(^2), (2.15) 

where the constant Scf? is 



_i f 2E 4^E ua 



(2.16) 

The new numerical constant is 



C4 ~ 0.282. (2.17) 

Wc note in passing that the logarithm in (2.16) represents the explicit subtraction of 
an effect analogous to a [quadratic x log] correction to the critical temperature for Bose 
condensation of dilute non-relativistic gases, discussed in Ref. [7]. In this analogy, the 
lattice spacing a in our simulation plays a role similar to the thermal wavelength in the Bose 
condensation problem, and u and Ji (which is chosen of order u) are proportional to the 
scattering length. 

In this paper, whenever we quote simulation results for a given value of ua and ra^, we 
are referring to simulations of the action (2.12) with parameters given by (2.13) and (2.15) 
with 0{- ■ ■) set to zero. When we quote values of A (0^)0, we will be quoting continuum 
values, as given by (2.15). This will be adequate to reduce the lattice spacing error to 0{a^) 
on individual measurements of A{(f)'^)c/u and to 0(a) on individual measurements Vc- 

The precise definition of continuum r, and its relationship to bare lattice ro, are in 
principle unnecessary if one's interest is only to determine A{(j)'^)c — one could simply find 
the critical value of Tq at any given lattice spacing and not worry about its relation to 
continuum definitions. However, as a practical matter, knowing the relationship facilitates 
using results at one lattice spacing to make a good initial guess of the critical value of Tq 
at a new lattice spacing. And, as discussed earlier, the continuum critical value of r is of 
interest in its own right. 

Throughout this paper, continuum r should be understood as defined by MS renormal- 
ization at a renormalization scale p,. That is, the continuum Lagrangian is the e ^ limit 
of the (3— e)-dimensional action 
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FIG. 1. The two fundamental UV-divergent graphs of continuum (j)^ theory. 



I' 



with the bare continuum parameter rbare related to the renormahzed r = Ty^{iJi) by 



and with 



'"bare — T -\- 



u 



(2.18) 



(2.19) 



(2.20) 



[The factor of e'^l"^ j^f^ in (2.20) is what distinguishes modified minimal subtraction (MS) 
from unmodified minimal subtraction (MS); the difference between the two schemes amounts 
to nothing more than a multiplicative redefinition of the renormalization scale. The constant 
7b = 0.5772 • • • is Euler's constant.] Three-dimensional 0^ theory is super-renormalizable 
and the only fundamental UV divergences of the continuum theory are those corresponding 
to the two diagrams of fig. 1. The first has a linear divergence, which is ignored by dimen- 
sional regularization. The second has a logarithmic divergence, which is the origin of the 
li^/e counter-term present in (2.19). 



C. Algorithm 

We evolve configurations in Monte Carlo time using a combination of site by site heatbath 
and multi-grid [8] over-relaxation updates. 

In section III, we discuss how we define a nominal value of the location Tc of the transition 
in finite volume. In order to scan for the transition, we need to be able to smoothly vary r. 
We use the standard technique of canonical reweighting. Having made a simulation at one 
value Tsim of r, and accumulated a Monte Carlo time series of values of 

0^^1|rf3^02(x), (2.21) 

expectations at nearby values r = rgim + can be calculated as 

exp (~\Vbr^\o\ 

V =xT^- (2.22) 

We estimate our statistical errors for the estimates of Tc and A(0^)c in a given simulation 
run using the jackknife method with 20 bins. The size of bins must be large compared to the 
decorrelation time, and this is verified in Appendix A, where we quote decorrelation times 
and sample sizes for each of our simulations. 
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III. FINDING THE CRITICAL POINT 



Systems only have sharply defined phase transitions in infinite volume, but we use the 
method of Binder cumulants [9] to obtain a good estimate, in finite volume simulations, of 
the critical value Tc of r. Specifically, we measure the cumulant 

as a function of r, where the 2-component vector is the volume average of (p: 

0=i jd?x<t^{^). (3.2) 

In infinite volume, the cumulant C is 1 in the ordered phase and 2 in the disordered phase. 
In large volume, there is a smooth transition between these two values, and the width of the 
transition region shrinks as the volume is increased. Specifically, for an L x L x L volume, 
the width in r of the transition region scales as L^^* in the L ^ oo limit, where yt = l/v 
and V is the correlation-length critical exponent. The value of the exponent is 

yt = Ijv ^ 1.49 (3.3) 

for the universality class of the 0(2) model. ^ 

One method of estimating the location of the transition in finite volume is to simply 
choose a fixed value C* of C between 1 and 2 and then define the nominal in finite volume 
as the r for which C(r) = C*. This leads to errors in (contributing to errors in other 
quantities measured at the transition) that would scale away as in the large L limit. 
This L^^* scaling of finite- size errors is typical of many prescriptions for defining in finite 
volume. 

The method of Binder cumulants improves the scaling of the finite-size error in from 
L~^* to L~^*~'^, where lo is the critical exponent associated with corrections to scaling, and 

ijj ~ 0.79 (3.4) 

for the universality class of the 0(2) model. ^ One version of the method is to measure the 
curves C(r) for two different large system sizes L\ and L2, and then estimate Tc as the point 
Tx where the curves intersect. Specifically, Binder [9] showed that, in the Li, L2 — > 00 hmit, 
the error scales as 



^ The value of roughly 1.5 comes from the scaling relationship v = {2 — a)/d and the fact that the specific 
heat critical exponent a is very small in this universality class. The best value of a is -0.01056(38) and comes 
from experiments in Earth orbit on superfluid He^ [10,11]. (See in particular endnote 15 of Ref. [11], and 
see also footnote 2 of Ref. [13].) For comparison, theoretical values (making use oia = 2 — dv as necessary) 
arc a = -0.011(4), -0.004(11), and -0.0146(8) from 3D series techniques [12], the e expansion [12], and 
Monte Carlo [13]. 



^Values are w = 0.789(11), 0.802(18), and 0.79(2) from 3D series techniques [12], the e expansion [12], and 
Monte Carlo [14]. 
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r^{Li, L2) - Tc 



1^2 — 



h-^ - 1 



(3.5) 



Lf - LI' 



1 - by* 



where b = L2/L1. Moreover, the value of the cumulant Cx at the intersection approaches a 
universal value Cc in this limit. ^ 

If one knew Cc in advance, then, for data in a given finite volume, a nice method for 
determining a nominal point of transition is to choose the r such that C(r) = Cc- The 
finite size error in Vc caused by this procedure also scales as L^^^^^ . This method is simpler 
and statistically a little cleaner than trying to find the intersection of two C(r) curves for 
two different values of L. This is the method we shall use, but first we need a value of Cc- 
Because Cc is universal, its value can be measured from simulations of any model in the 
same universality class. We will use a value determined by Campostrini et aL [13], who also 
simulate a 2-component lattice (p"^ theory. They obtained Cc = 1.2430(1)(5), where the two 
numbers in parenthesis represent their statistical and systematic errors, respectively. In the 
rest of this paper, when quoting results for the transition in a given finite volume, we will 
mean the point where 



for that volume. As a check, we have also made a much less accurate determination of Cc 
using our own simulations, which is discussed in Appendix D. We find Cc = 1.2402(7), 
which we suspect has a systematic bias, discussed in the Appendix. The difference between 
these values of Cc is insignificant for our purposes. We have checked that the difference 
would only change our final results by a tiny fraction of an error bar. 



As discussed in the introduction, the only parameter of the continuum problem at the 
critical point is u- So the only length scale of the problem is 1/u- The relevant measure of 
the size L of a finite- volume lattice relative to this scale is therefore Lu (and similarly the 
relevant measure of lattice spacing is ua). Fig. 2 shows a plot of our results for A(0^)c on 
L X L X L lattices vs. Lu for ua = 6. As can be seen, our largest Lu values are clearly large: 
the data clearly shows nice convergence towards an infinite volume limit. To understand 
the size of the remaining finite- volume error, we will want to fit the volume dependence at 
large Lu to an appropriate scaling law. 



The scaling of large L corrections depends on universal critical exponents of the 0(2) 
model, which is in the same universality class as the classical N = 2 Heisenberg ferromagnet. 
The language of critical exponents in the 0(2) model is borrowed from correspondence with 



^For a nice numerical demonstration of this universality, see the Monte Carlo studies of Ising universality 
in Ref. [15]. 



C(r) = 1.243 



(3.6) 



IV. VOLUME DEPENDENCE OF A{(f)^)c 



A. Large volume scaling of A(0^)c 
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FIG. 2. Simulation results for A{<lP')c/u vs. lattice size in physical units (Lu) for ua = 6. The large Lu 
curve shown is a fit to the rightmost 4 data points, with confidence level 61%. 
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the Hciscnbcrg magnet, and t = r — Tc is referred to as the reduced "temperature" in this 
context, C oc {In Z )/ dt^ as the "specific heat," and so forth. It is important to emphasize 
that this language holds the potential for confusion because, in many applications (such as 
Bose-Einstein condensation at fixed density), t = r — Vc is not, in fact, finearly related to 
the actual temperature of the system near the critical point, and d'^{\nZ)/dt^ is not directly 
the physical specific heat.^ In any case, the analog of the "energy density" E/V under this 
correspondence is 

where V — is the system volume and d — 3 is the dimension. The problem of under- 
standing how the finite volume corrections to (0^) scale with L can therefore be stated in 
this language as the problem of how corrections to the "energy density" scale with L for 
anything in this universality class. 

The standard finite-size scaling ansatz provided by renormalization group methods is 
that the free energy density / = —V"^ In Z scales as^ 

fit, {uj}, L-') = Aeg(t, {uj}) + h-<'U^^{hy% {u^hy^}, b/L) (4.2) 

for periodic boundary conditions, where /reg and /sing generate the so-called regular and 
singular parts of the free energy in the infinite volume limit. (See Refs. [17,18] and references 
therein.) The length b is an arbitrary renormalization scale (block size), and {uj} denotes 
the set of infrared-irrelevant operators (with corresponding yj < 0). Standard notation 
for critical exponents is u = l/yt, and we will denote the smallest \yj\ associated with the 
irrelevant operators {uj} as cu. 
Choosing b — L, 

fit, Uj, L-') = f,,^{t, Uj) + L~''f,^^{Ly% Ujiy^, 1). (4.3) 

The usual infinite- volume scaling form is obtained by taking L — > oo with t fixed, which, for 
the limit to exist, requires 

/sing(T,0,l)~T'^/^* as T^OO. (4.4) 

In contrast, for fixed L, the free energy will be analytic in t, since there are no phase 
transitions in finite volume. We can make a Taylor expansion of the finite- volume free energy 



uniform, non-relativistic Bose gas is a constrained system: the particle density n is fixed. This constraint 
causes a different relationship of model parameters and the physical temperature than for unconstrained 
systems [16]. For example, the critical exponents x = (a, /3, 7, />) of the actual system are related to the 
standard exponents x = (a./?. 7,;/) of the field theory by (i) a = — — a), and x = x/{l — a) for the 
others, if a > 0, or (ii) i = a; if a < 0. This relation explains the difference between mean-field theory 
exponents for the 0(2) model (e.g. a = 1/2) and the exponents of a non- interacting Bose gas (e.g. a = — 1). 

^We have not included a "magnetic field" h coupling to (p, with a corresponding argument b^f^h in /sing, 
because we will only be interested in the case h = and are not interested in correlations of cj) (as opposed 
to 0^), which could be generated by derivatives with respect to h. 
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(4.3) in t (as well as {ujL^^}), which should be a Tiseful expansion when the arguments of 
the scaling piece /sing are small. That is, for situations where L^H ^ as we take L ^ oo, 
we can Taylor expand (4.3) as 

/ = (Ao + BoL-" + CoL-'^-) + • • •) 

+t(Ai + BiL-'^+y' + CiL-'^+f*-'^ + • • •) 
+t\A2 + B2oL-''+^y* + ---) 

+ ■■■, (4.5) 

where we have only displayed the leading corrections due to irrelevant operators. Differen- 
tiating with respect to t to get the energy density, we find that (0^) scales in large volume 
as 

(02^ = (Ai + BiL-'^+y' + CiL-'^+y'-^ + ■■■) 

+2t{A2 + B^oL-'^+^y^ + ...) + .... (4.6) 

As mentioned earlier, use of the method of Binder cumulants to determine means that, 
in our application, 

t ~ L-^*-'". (4.7) 

This indeed satisfies the condition L^H ^ as L — >■ oo, and so the expansion (4.6) is 
appropriate. For our application, we then have 

(02^ = + BiL-'^+y* + A'^L-y*-"" + c'lL-'^+y*-^ + ■■■. (4.8) 

Using the standard scaling relation a — 2 — vd for the specific heat scaling exponent a, this 
can be rewritten in the form 

{(f)^) ^Ai + B^L-^^-°'^y' + A^L-y'-"" + c;l-(^-")2'*-'^ + • • • . (4.9) 

The value of a in the three-dimensional 0(2) model is very small: a ~ —0.013, corresponding 
to the value yt ~ 1.49 quoted earlier [10-12,14]. 

The renormalizations Z,f. and 5(jP' that convert (0^)iat into A{(f)'^) in (2.15) do not introduce 
any new powers of L, and so the form of the large-L expansion of A(02) is the same as that 
for {(p'^) in (4.9) above, though the coefficients are different. 



B. Large volume scaling for a — 

Because we do not have large volume data spanning many decades in L, a is zero for 
all practical purposes. And so one might as well use the a = limit of the large-volume 
scaling (4.9), which corresponds to yt = d/2. Typically, a = generates logarithms in an 
RG analysis, which can appear as a superposition 

lim = qs^ In s (4-10) 

a-»0 a 
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of power laws 5^+"?" and for some variable of interest s. One might therefore expect that 
the a ^ limit of the large- volume scahng (4.9) is 

(02) ^^i + SiL-'^/^ + L-'^/2-'"(C In L + £>) + •••. (4.11) 

It's useful to verify the presence of a logarithm. If it weren't there, we could include the first 
corrections to scahng in large L fits of our data using a 3-parameter fit [Ai, A2, D] rather 
than a 4-parameter fit [^li, Bi, C, D; or equivalently Ai, Bi, A2, C[ in (4.9)]. 

The existence of the logarithm can be directly related to the well known logarithmic 
divergence of specific heat with t when a = 0. We give a renormalization group analysis in 
Appendix C that makes this explicit. Here, let us just note that the logarithm follows from 
an old proposal by Privman and Rudnick [19]. Ignore corrections to scaling for the moment, 
and suppose that at o; = the free energy had the general form (4.3) discussed earher: 

fit, L-') ~ /reg(t) + L-'U.itL"/'). (4.12) 

In order to get a logarithmic divergence ln(i~^) of the specific heat in the infinite volume 
limit, we need a term ln(t~^) in the free energy in that limit. So we must have 

/sing(T) ~ ATMn(T-^) as T^OO. (4.13) 

But this would give / ~ /reg(^) + ln(t^^L~'^/2), which doesn't have a good L ^ 00 limit. 
The solution is to suppose that the a — version of the free energy is instead 

fit, L-') /r,g(i) + At' IniL'/') + L-Vsing(iL<^/'), (4.14) 

which is what Privman and Rudnick proposed. Note that the new term is analytic as t ^ 
for L fixed, as it must be. The new term gives a tlnL contribution to the energy density 
E/V , which, for the t ~ L-y^-^ of interest to us, gives rise to the logarithm term in (4.11). 



C. How large is large volume? 

Before proceeding to numerical fits of the large volume dependence, it is useful to first 
have an idea of how large L should be before one might reasonably hope for large volume 
scaling to hold. From Fig. 2, we see that the finite volume corrections to the continuum 
value of A(02)c become 100% where the data crosses zero, at roughly Lu ~ 200. So one 
might guess that this is very roughly the scale where scaling starts to set in — the scale that 
separates short-distance perturbative physics from long-distance critical-scaling physics. We 
can check this rough assessment from the other side. In the hmit of small Lu, the physics 
of fluctuations is perturbative, and one can analytically compute the expected value of 
A(02)c order by order in powers of Lu. We perform this computation in Appendix E in the 
continuum (wa = 0) to next-to-leading order (NLO) in Lu, with the result that 

A(02^c _ 6.44003 0.451570 , r^nT. .\-i/2 



u 



(Lm)3/2 Lu 



+ 0{iLu)-^") (4.15) 



for our critical value (3.6) of the Binder cumulant. The resulting curve is shown in both 
figs. 2 and 3 (as well as the leading-order small Lu result in fig. 3). One can see that the 
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Lu 

FIG. 3. Simulation results for A{(p'^)c/u vs. lattice size, showing more small volume data than Fig. 2. 
The dot-dash curve shows the next-to-leading order small Lu expansion of (4.15). The dotted curve shows 
just the leading-order term of that expansion. 

small Lu expansion becomes unreliable past Lu ~ 200, which is the same as the previous 
scale estimate. (There is obviously nothing precise about this statement. The scale at which 
the second term in the small Lu expansion becomes 50% of the first term, for example, is 
Lu ~ 50.) 

Interestingly, the scale Lu ~ 200 is close to what one might estimate on the back of 
an envelope from a large-// approximation to the theory. In large N, one replaces the 
0(2) theory studied here by an 0(iV) theory of scalar fields and solves the theory in 
the approximation that N is large — a program pursued for the problem of Bose-Einstein 
condensation of a non-relativistic gas in Refs. [20,21]. In the N ^ oo limit, one introduces 
an auxiliary field a whose propagator represents a geometric sum of bubble diagrams, such 
as shown in fig. 4. The corresponding resummed propagator is proportional to 

1 

where 5]o(p) represents the basic massless bubble integral 

The propagator (4.16) changes from its large momentum behavior (^ constant) to its small 
momentum behavior (oc p) at a scale given by 3/Nu ~ So(p). This corresponds to p ~ 
Nu/AS and so distance scales of order L ~ 27r/p ~ ddn/Nu. Setting N — 2, we find 
Lu ~ 487r ~ 150. Of course, one would never hope that an estimation this crude would be 
useful beyond, at best, the factor of 2 level. 
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(4.16) 



(4.17) 



>-< = xoo<^>oo< 



^ >ooo< ^ - 

FIG. 4. Bubble chains, which are the source of scale dependence in large N calculations of A{(p)c in 
0{N) theory [20,21]. 

Our discussion of system size has implications for the vahdity of an early numerical study 
of the critical temperature for Bose- Einstein condensation for non-relativistic gases [22]. We 
discuss this in Appendix F. 



V. NUMERICAL EXTRAPOLATION OF Lu 



oo 



A. Extrapolating A(</)^)c to Lu oo 

We now examine fits of the ua=6 data of Fig. 2 to the large L scaling form of (4.11). The 
circles in Fig. 5 show the result of extrapolating to L — cxd using the leading-order scaling 
corrections in (4.11) — that is, ignoring the sub-leading terms, which involve the exponent cu. 
For comparison, the diamonds show what happens if we ignore finite-size effects altogether; 
the corresponding confidence levels are terrible. The vertical dashed line in the figure marks 
the moderate system size L = 200 (discussed in the previous section), below which you 
should become suspicious of any attempt to fit the data to a large-volume scaling form. 

A simple method for assigning a final result for the extrapolation is to take our best fit 
with a correct scaling form. As seen in Fig. 5, the leading-order scahng form A + BL~'^/'^ is 
perfectly adequate for fitting our large Lu data. Our procedure is to fit the data for lattice 
sizes greater than or equal to a given Lmin, decreasing this minimum size for as long as 
the fit remains stable with a reasonable confidence level. We take as our estimate the 61% 
confidence level fit to Lu > 384, which gives 



U 



-0.001289(9), 



(5.1) 



«a=6 



which is depicted by the shaded region of Fig. 5. 

As described in Rcf. [1], we will actually use Lu = 576 as a reference point from which 
to derive finite volume and finite lattice spacing corrections. Fig. 6 is similar to Fig. 5 but 
shows the size of the finite-size correction at Lu = 576, as determined by the fit. The best 
fit (the 61% confidence level one) gives 



rA(02)ci 






u 


Lu=576 


u 



= 0.000241(7), 



(5.2) 



Lu— >oo 



and we will use this difference, rather than the limit (5.1), in what follows. [The difference 
is determined by the coefficient B of our fit to -|- BL~'^/'^^ and so is determined by all the 
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FIG. 5. Extrapolations of the ua=Q data of Fig. 2 to X — > oo under various assumptions of the functional 
form. The horizontal axis shows the smallest value of Lu of the data used for the fit. Percentage confidence 

levels are written next to each extrapolation whenever non- negligible (> 1%), and "NA" (not applicable) 
refers to those cases where the number of data points used for the fit equals the number of parameters. 



data points of that fit; it is not simply our {Lu, ua) — (576, 6) data point minus the hmit 
(5.1), which would produce a larger error.] As wc shall see, wa = 6 is a reasonably small 
value of ua, and we expect this to be a good estimate to the finite volume corrections in the 
continuum {ua — > 0) limit. 

As a check that corrections to scaling will not radically alter our results, we show as 
squares in Fig. 7 the result of fits to (4.11). The values are consistent with the previous 
result, with larger errors because we are fitting more parameters. The triangles show a 
simplified fit, with one less parameter, that ignores the logarithmic dependence. 



B. Extrapolating Tc to Lu ^ oo 

We will now make a similar analysis of finite-size effects for the critical value Vc of r. The 
continuum value of r is convention dependent — it depends on one's choice of renormalization 
scheme and renormalization scale. As discussed earlier, our convention will be to define r 
with MS regularization at the renormalization scale fi = u/3. The conversion formula (1.6) 
to other choices of p, can be extracted from (2.13c) and the fact that the theory is super- 
r enormalizable . 

Fig. 8 shows, for ua — 6, the dependence on system size of our finite- volume determina- 
tions of rc{u/2,)/u^. As discussed in section III, the finite-size corrections are expected to 

scale as L~^^~'^ as L ^ oo. Fig. 9 shows the result of extrapolating an infinite- volume result 
by fitting various subsets of the data to A-\- BL~'^'^/'^^~'^ . Taking the highest confidence level 
fit. 
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FIG. 6. As Fig. 5 but shows the magnitude of the finite-size correction to A(0^)c/u at Lu = 576, as 
determined by the fit. 



T3 
C 
(0 
CO 

LO 



o 

CD 



CD 



CD 
O 

c 

CD 
i_ 

CD 



0.0005 
0.0004 
0.0003 
0.0002 
0.0001 


-0.0001 



%4 -- 

S q 

6 31 3^ 30 
46 33 49 



□ a+bL ^'^+cL ^'^ "(In/.+d) 

. 1 -1.5 1 -1.5-(0 

^ a+bL +cL 



NA 
36 



NA 



500 

smallest Lu used for fit 



co=0.79 

1000 



FIG. 7. As Fig. 6 but showing extrapolations that include corrections to scaling. Where two confidence 
levels are listed on top of each other, the upper one is for the square data point and the lower for the 
triangular data point. 
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FIG. 8. ua = 6 data for the nominal value of r^jvi} (defined in MS renormalization at renormalization 
scale /z = u/3) as a function of system size. The fit shown by the dashed line is the large-volume fit to all 
but the leftmost data point. 



rc(M/3) 



0.0028828(6). 



(5.3) 



- ua=6 



Fig. 10 shows the size of the finite volume error at our canonical Lu — 576, as determined 
by the fits, 



rc(w/3) 




rc(w/3) 






Lw=576 




Lu=oo 



= -0.00000604(26). 



(5.4) 



VI. NUMERICAL EXTRAPOLATION OF ua ^ 

A. Extrapolating A(^^)c to ua ^ 

We will now discuss numerical results for the continuum limit a — > while holding the 
physical volume of the lattice fixed. This can be expressed as — > holding Lu fixed. 
Since the number {L/cbf of lattice sites we can practically include in a simulation is limited, 
we can obviously explore smaller values of ua when we fix smaller values of Lu. As a test 
that we understand our lattice spacing errors, we have therefore made several simulations at 
the rather moderate system size of Lu = 144 (see the discussion of sec. IV C). The results 
for the dependence of A(0^)c/m on ua are shown by the squares in Fig. 11. If our corrections 
for 0{a) errors have been calculated correctly, the remaining error should be O(a^). Indeed, 
all but the largest two ua data points in Fig. 11 fit very well the functional form A + B{uay, 
with confidence level 94%. 
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FIG. 9. Extrapolations of the ua=6 data of Fig. 8 to L ^ oo using the functional form A + 
The horizontal axis shows the smallest value of Lu of the data used for the fit. Confidence levels are written 
as described for Fig. 5. 
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FIG. 10. As fig. 9 but showing the magnitude of the finite-size correction to rc/u^ at Lu = 576. 
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FIG. 11. The squares show resuhs for A((/)^)c/u vs. ua at Lu = 144. The hnc through them is a fit of 
the first 6 points to A + B{ua)^ . The diamonds represent the corresponding uncorrected data, as described 
in the text. To guide the eye, a straight hne has been fit through the first three points, though the actual 
smah ua behavior is hnearxlog. 

It's interesting to compare to what would have been obtained if we had used the same 
data but instead plotted the uncorrected 



(6.1) 



27ra. 

vs. Uoa, where $ and Uq are the bare lattice fields and coupling of (2.1). This uncorrected 
data is represented by the diamonds in Fig. 11. One clearly sees the 0{a) corrections. We 
should make clear that this is still Lu = 144 data and is not Luq = 144 data, which would 
have required additional simulations. 

Fig. 12 shows the ua dependence at a reasonably large physical system size of Lu = 576. 
We show extrapolations of the ua ^ limit in fig. 13. Here, the 10% confidence levels of 
fits to A + B{uay are less spectacular than the Lu = 144 data, though not unreasonable. 
Because of the lower confidence levels, we have been a little more conservative in our error 
estimate. We take as our result for the ua ^ limit the shaded region of fig. 13, which has 
been chosen to cover both the 10-15% confidence level fits: 



u"^ 



= -0.000957(15). (6.2) 

Lu=576 



Combining this with the finite volume correction (5.2), and adding errors in quadrature, 
we arrive at our final value (1.4) for the infinite-volume continuum limit. Note that the 
dominant error in this estimate comes from the ua —>■ extrapolation. 

B. Extrapolating Vc to ua —> 

The circles in Fig. 14 show the dependence of Vc on ua at the medium system size of 
Lu = 144, where we can simulate down to relatively small values of ua. The data fits well 
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FIG. 12. Results for A(0^)c/m vs. ua at Lu = 576. The line is a fit of A + B{ua)^ to all but the rightmost 
data point. 
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FIG. 13. Extrapolations of the data of Fig. 12 to ua = fitting to the form A + B{ua)^. The horizontal 
axis shows the maximum value of ua used for each fit. Confidence levels are written as described for Fig. 5. 
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FIG. 14. Results for Tc/v? vs. ua at Lu = 144 and Lu = 576 (defined in MS renormalization at 
renormalization scale p, = u/3). The lines are linear fits. 

to a linear dependence on ua, with the fit shown to the first six points in the figure having 
a 38% confidence level. The triangles in the same figure show similar dependence for the 
reasonably large system size of Lu — 576. The corresponding extrapolations to ua — based 
on linear fits are shown in fig. 15. We take the ua — extrapolation to be 



rc{u/3) 



u^ 



= 0.0019141(21), (6.3) 

L«=576 



as shown by the shaded region of the figure. Combining with our estimate (5.4) of the finite 
size error at Lu = 576 gives our final result (1.5) for the infinite- volume continuum limit. 
Again, the dominant error comes from the ua ^ extrapolation. 



C. Interdependence of Lu oo and ua extrapolations 

In this paper, we have treated the extrapolation of our finite vohimc correction, taken 
from ua=6 data, as independent from our continuum extrapolation, taken from Lu=576 
data. Consider the case of A(0^)c/m. In principle, the two extrapolations are not com- 
pletely independent because the coefficient B in the form A + BL~'^/'^ used for our large 
L extrapolation is not universal. Its value will have some dependence on ua. For our pro- 
cedure, we would like to know the value of B oX, ua = 0, to obtain the large L correction 
(5.2) at ua = 0. Instead, we have the large L correction at ua = 6, which will be slightly 
different. In numerical terms, we will have underestimated our systematic error if B[ua=6] 
differs from B[ua—0\ by ~ 7%, which is the final error on A(^^)c/u quoted in (1.4) relative 
to the size of the finite volume correction (5.2). 

How does the coefficient B depend on ua as ua ^ 0? Because of our improvements to 
the action and to the operator 0^, the answer is that small ua corrections to B vanish as 
{ua)^. To see this, note that the only quantities we have not matched through 0{{uay) are 
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FIG. 15. Linear extrapolations of the Lu = 576 data of Fig. 14 to ua = 0. The horizontal axis shows 
the maximum value of ua used for each fit. Confidence levels are written as described for Fig. 5. 

r and 50^. However, neither of these contribute to B. Scf)^ is a constant and therefore Lu 
independent, and our procedure tunes to the critical value of (the continuum parameter) 
r whether or not the lattice parameter is correctly matched to the continuum one. The 
residual difference between lattice and continuum r does appear in Eq. (2.16), but, because 
r — Tc vanishes as jj-'^-'^l'^ ^ this does not contribute to B either. We should not expect such 
a high order correction to be of any consequence. 

There is an alternate way of fitting our data which makes it easier to check this. We 
can fit all data at large Lu and relatively small ua to a form which allows both spacing and 
volume dependence, 

^^ = A + S(LM)-2'* + C(«a)2. (6.4) 

Based on our fits at fixed Lu and fixed wa, we expect this form is sufficient if we use only 
data with Lu > 384 and ua < 8. Fitting all such data gives A{(l)^)c/u = -0.001194(8) with 
X^/d.o.f. = 15.5/8 (5% confidence level). The result is consistent with our quoted result, 
and has smaller errors. The quality of the fit is not very good, which is why we quote the 
larger error bars of our other analysis. However, adding a term which allows B to vary with 
ua (that is, a ua and Lu dependent term) does not improve the quality of the fit. Adding a 
new term D{Lu)~y*{ua)"^ and doing a 4 parameter fit changes by less than 1 for m = 2 
or m = 3; so the data show no evidence that such a term is large enough to be important. 

We do not understand why the fit to all data has such a large (x^/d-O.f = 15.3/8). 
The three lowest ua data points are the main outliers; the two ua = 3 points are each at 
2a. However, the {Lu, ua) = (384, 3) point is off in the opposite direction as the {Lu, ua) = 
(576, 3) and {Lu, ua) = (384, 4) points, so there is no systematic trend in the residuals and 
we have no reason to believe that some additional lattice spacing or volume dependent term 
is missing in our fitting procedure. 

Finally, there is a minor, independent issue in our main analysis. The statistical part 
of our errors in our original infinite volume and continuum extrapolations are correlated. 
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since both of these extrapolations included the {Lu,ua) = (576,6) data point. However, 
since the error from the continuum extrapolation dominated that from the infinite volume 
extrapolation, any statistical cross-correlation between those errors will not have a significant 
effect. 
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APPENDIX A: TABULATED DATA 
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To give a rough idea of the size of our data sets, we have listed above a nominal decorre- 
lation time r for each simulation, along with the amount of data we have in units of 2t. In 
our convention, r = 0.5 represents completely uncorrelated data. Our decorrelation times 
arc in units of sweeps, where one "sweep" consists of both a heatbath sweep and a multi- 
grid update. More than one time is shown in cases where simulations were made at different 
values of r (before reweighting). 

Our nominal decorrelation time is the largest integrated decorrelation time of the various 
expectations required in the computation of the Binder cumulant and of (0^) . These include 
the integrated decorrelation times associated with measurements of P, -P0^, and 
(and also times associated with cross-correlations between these), where 
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P = exp(-|y5r(/)2) (Al) 

is the canonical reweighting factor, as in (2.22). We found in all cases that the longest 
decorrelation time was that associated with 0^. The integrated decorrelation time for a 
single operator is given by 

1 - CH 

" 2 ^ ntl C{0) 

where 

1 N~n / 1 ^ \ ^ 

C{n) = g - g O.j , (A3) 

is the auto-correlation function associated with the desired operator O. In practice, the sum 
in (A2) must be cut off because of degrading statistics, and we cut it off when C(n)/C(0) 
first drops below 0.05. 



(A2) 



APPENDIX B: MATCHING CONTINUUM AND LATTICE THEORIES 

1. General discussion 

In this appendix we discuss the 0{a) and O(a^) improvement of three-dimensional scalar 
field theory on the lattice. For the sake of generality, and because it is not any harder, we 
will discuss O(A^) scalar field theory of real fields = (0i, 02, ■ • ■ (/)n)- The case of interest 
to the present work is = 2. As in (2.12), the lattice Lagrangian is defined to be 



X 

where 



(Bl) 



N 



<I>'^E<I>1 (B2) 



and (j), r, and u are just the (UV renormalized) continuum fields and parameters in lattice 
units. So 0iat = a^^Vcont, «iat = a«cont, and riat = aVcont- Of these, rcont is the only contin- 
uum parameter that requires UV renormalization and should be understood as renormalized 
with dimensional regularization and the MS scheme. In this appendix, we will calculate, 
to a given order in lattice spacing, the necessary counter-terms Z^, Zr, 5r, and 5u required 
to implement this correspondence between continuum and lattice variables. To match the 
lattice and continuum theories to high orders in a, we would need to include other operators 
in our lattice theory, such as 0^, 0^|V0p, and so forth. However, these will turn out to be 
unnecessary at the order to which we will work. 

A given local lattice action will never perfectly reproduce the continuum action. For 
small ua, the error in how a given lattice action treats physics at the distance scale a can 
be computed and compensated for by perturbative calculations. The discrepancy in how 
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(c) (d) 

FIG. 16. The one loop graphs needed for the renormahzation. A eross represents a cfP" insertion, which 
in turn could represent either (i) a perturbative insertion of the rtj)^ term of the action or (ii) the operator 
associated with calculating the expectation {(f^). 

a given lattice action, with given parameters, treats the non-perturbative physics at the 
distance scale 1/u cannot. It is therefore important that the lattice action be close enough 
to the continuum action that errors at the scale 1/u are higher order in a than whatever is 
desired. In order to improve our simulations and measurement of A(0^)c to 0(a) accuracy, 
it is necessary to (a) match the lattice action and parameters so that, at the scale it 
reproduces the physics of the continuum action up to and including 0(a); and (b) match the 
lattice and continuum definitions of operator cfp' to the same order, so that the measurements 
of A{(f)'^)c will match up. Actually, the first requirement is slightly overstated. To measure 
A(0^)c to 0(a), it is not necessary to match the lattice and continuum parameters r to that 
order. That's because, in the simulations, we will simply vary the coefficient of 0^ in the 
action until we find the transition - we don't need to know its relation to rcont to do this. 

In this appendix, we will match rcont to the lattice just to 0(a°). This will make possible a 
determination of the continuum value of the critical value Vc, but 0(a) errors will remain and 
must be removed by extrapolation to the continuum limit. We will match the continuum 
and lattice definitions of the operator 0^ through 0(a), so that we can make an 0(a) 
improved calculation of (0^)c- We will also match the lattice action through O(a^) [except 
for the matching of r just discussed], rather than simply 0(a), because it is not that difficult 
[compared to the 0(a) matching of 0^] and should improve the ua dependence of infinite- 
volume extrapolations of our data, as discussed in section VI C. 

At tree level, the infrared behavior of the lattice and continuum theories are the same 
up to corrections suppressed by at least a^; the power of a can be higher if we choose the 
lattice appropriately. However there are "radiative" corrections induced in the IR physics 
by the nonlinearity of the theory, together with the UV difference between the lattice and 
continuum theories. The diagrams of Fig. 16, for instance, differ on the lattice from their 
continuum values because the dispersion relations of the scalar fields differ in the UV, and 
the loop momentum integration integrates over this region. These effects are suppressed by 
powers of the couplings; the radiative correction to the scalar 4 point function, for instance, 
clearly depends on m^. On dimensional grounds, the difference between lattice and continuum 
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FIG. 17. A diagram contributing to the mixing of the cjp' and <j)^ operators. 

values of this diagram must go as w^a, so these diagrams lead to 0{a) differences between 
the lattice and continuum theories. The difference between lattice and continuum values of 
these diagrams can be removed by a renormalization of the parameters of the lattice theory. 

As discussed in the main text, we will be interested not only in how to renormalize the 
parameters of the action but also in how to translate expectation values of the operator (f? 
between the lattice and continuum. It will be convenient to talk directly about the operator 
0^, which is associated with UV divergences in the continuum, rather than A(0^), which 
is not. We will define 0^ in the continuum also using MS renormalization. In general, 
operators with the same symmetry can mix under renormalization, and the lattice operator 
(fp' will correspond to some superposition of the unit continuum operator, the renormalized 
0^ continuum operator, and higher-dimensional renormalized continuum operators such as 

a~^(0^)lat = Co + C2(0^)cont + C4(0^)cont H ■ (B3) 

However, our particular interest is in expectation values at the transition. For this applica- 
tion, the effects of higher-and-higher dimensional operators on the right-hand side of (B3) 
become suppressed by more and more powers of ua. For example, C2 is 0(1), and the expec- 
tation value of the renormalized continuum operator cfp' is, by dimensional analysis, 0{u) at 
the transition. So the C2(0^)cont term contributes 0{u) to the expectation. In contrast, the 
lowest-order diagram that contributes to mixing between 0^ and 0"^ is shown in Fig. 17 and 
gives C4 = O(M^a^), where the v? counts the two vertices in Fig. 17 and the a? then follows 
by dimensional analysis. But, also by dimensional analysis, the renormalized continuum ex- 
pectation (0^) is Oiv?). So the C4(0^)cont term contributes 0{u^a^) to the expectation, down 
by three powers of ua from C2(0^)cont- Our goal will be to compute the 0{a) corrections to 
A(0^)c, for which it is therefore adequate to compute just Cq and C2 above. 

The diagrams which contribute to C2 are the same as those which contribute to the 
multiplicative renormalization of the coefficient r of 0^ in the action, and l/c2 is the same 
as the Zj. introduced previously. Rearranging the terms in (B3), we will write 

a(0')cont = Zr{4>%.i - 50' + 0{{uaf), (B4) 

where represents a c- number (aco/c2) parameterizing mixing with the unit operator, 
which we shall calculate. 
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We will see later that, for the purpose of calculating the 0{a) corrections to A(0^)c, all 
that is strictly required is one-loop results for Z^. Z.,.. 6u. and 5r, and three-loop results for 
To determine Tc [just at 0(a°)] requires tree results for Z^, and Z^, and two-loop 
results for 5r. If one made a three- loop computation of 5r (which we have not), one could 
obtain the 0{a) corrections for Tc as well. To match the action through O(a^), except for r, 
one wants a two-loop determination of Z^ and 8u. As discussed in section VI C, it will also 
be useful to have a two-loop result for Z^. 

Now we turn to the calculations. For the remainder of this appendix, we will work 
exclusively in lattice units (a = 1). 



2. One-loop results 

A one loop renormalization calculation will determine the 0(a) contributions to the 
quantities Z<^, and Z,., and will find the 0(l/a) contributions to and The 
details of the power counting used here can be found in [23]. In the small lattice spacing 
limit, the momentum scale r^^^ associated with the parameter r is small compared to the 
scale 1/a where the lattice and continuum theories differ. So, for the specific purpose of a 
calculation to match the lattice and continuum theories, the term in the action (as well 
as the term) may be treated perturbatively. 

The required graphs are shown in Figure 16. Evaluating the graphs requires choosing 
a lattice Laplacian, and for completeness we will consider both the unimproved (2.2) and 
improved (2.3) choices described in the main text. The evaluation of the one loop graphs 
requires two integrals: 

^= / ^1 fB5) 

47r ~ Jbz (27r)3 fc2 ' ^ > 

[ ^ J__ / ^1 fB6l 
47r kz (27r)3 (fc2)2 (2T:f ' ^ ^ 

Here we use the shorthand BZ to mean k lies within the Brillouin zone, meaning each 
ki e [— 7r,7r]. The notation k'^ is introduced in Eq. (2.4). The integrals which determine ^ 
are each IR singular and some regulation is implied, for instance adding to both k^ and 
/c^ and taking the limit as — > 0. The numerical values of the integrals, accurate to ±1 in 
the last digit, are ^ 

Eu = 3.17591153562522 , Ei = 2.75238391130752 , 

= 0.152859324966101 , 6 = -0.083647053040968 . (B7) 

Note that the sign of ^ depends on whether we use an improved lattice Laplacian. This 
is possible because ^ represents the difference of a graph between lattice and continuum 
theories. The lattice contribution is larger inside the Brillouin zone, but the continuum 



^An analytic result [6] for Sy is f (18 + 12^2 - 10\/3 - 7\/6) [K((2 - \/3)2(\/3 - \/2)2)] ^, where K is the 
complete elliptic integral of the first kind. 
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FIG. 18. All required two loop grapliH and one loop graphs with one loop counterterm insertions, shown 
as heavy dots on lines (mass counterterms) or at vertices (coupling counterterms). The last eight graphs 
cancel in pairs. Diagrams (a), (d), and (e) are not separately IR convergent; diagram (d) must be distributed 
between the other two to produce IR convergent integrals. 



integral receives contributions from outside the zone as well; the sign depends on which 
effect is larger. 

At one loop the renormalizations are (in lattice units) 
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3. Two-loop results 



For the renormalizations and Z^, it makes no sense to carry the renormalization 
to two loops unless we use the improved lattice Laplacian, because already at tree level 
the unimproved Laplacian gives O(a^) level errors in the propagator. The two loop results 
require several more graphs as well as the inclusion in one-loop graphs of one-loop mass 
and coupling counterterms. (See Figure 18.) Four more integrals are needed, which we will 
evaluate in a moment. The complete 2- loop result for general A^" is 
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The three numerical constants Ci, C2, and C3 are given for the improved Laplacian in (2.14). 
Only C3 can be usefully defined for the unimproved case/ where it is Cs^u = .08848010. 

Now we detail the calculation of the constants Ci through C3. We will use the following 
shorthands: /^ bz means / (i^/c/(27r)^, with range the Brillouin zone [— tt, tt]^; whereas /j^. jj3 
is the same but integrated over all 3-space. Further, we define the following integrals, which 
will come up repeatedly: 



h{p) 

Hp) 



q,BZ q'^{p + q)^' 
1 



(B18) 
(B19) 



We begin with the two-loop vertex correction, graph (a) of Figure 18. The required 
integral, including the appropriate amount of the one loop counterterm graph (d), is 
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(47r)2 Jk,BZ (A;2)2 
We re-arrange the original integral into three parts. 
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Ic{k). (B20) 



(B21) 

fc,5R3_BZ 8k^ ^ ' 




All three of the above integrals are convergent, provided we use the improved Laplacian. 
The first integral is IR well behaved because the two counterterms cancel li,{k) up to a 
correction, which in the small k limit is K{k) - (1/8/c) - (C/47r) ^ .0125438 /cV^tt. 

To get accurate numerical answers, we perform all 3-D integrals by quadratures. Dealing 
with the double poles which appear in /l is touchy, and requires adaptive mesh refinement 
techniques. We improve the precision of each final result by repeating the full integration at 
several spacings and extrapolating (Richardson extrapolation). The first integral in (B21) 
gives .0360003/1671^, and the second gives .054568958/1677^. The last integral, over 3fJ=^-BZ, 
can be re-arranged into 
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(B22) 



''In ref. [24], this constant is called C- 
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which can be performed very accurately and gives —.035507296027 • • • /IGtt^. These sum to 
give Ci = .0550612. 

Besides this graph, there is graph (e), which gives 
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1 



which is not IR convergent. However, including -2 times the counterterm diagram (d). 



(B23) 



k,BZ {k'^Y ) \Jk,BZ {k^) 



(B24) 



gives — (^/47r)^. No new integrals are required. It is a nontrivial check on the calculation 
that the sum of the coefficients arising from diagrams (a) and (e) precisely absorb diagram 
(d). 

The next integral is the O(p^) contribution from the setting sun diagram (b). 



The first trick is to note that 



Jk,BZ \ {k+pf k'^J 



(B26) 



just by shifting the integration variable for the first term. So we may add an arbitrary con- 
stant to Ii{k) in (B25), and we choose the constant — ,^/47r. This will prevent IR divergences 
in what follows. We are now free to expand l/{k+pY to second order in p. After averaging 
over directions for p, we find 



{k+pf 



]_ / 8 sin fcj— sin2fcA ^ i Z' 
3 3 J 3 \ 



4 cos fcj— cos 2k 
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p^M{k). (B27) 



The equivalent expression in the continuum case is p'^/Sk'^. Re-arranging the terms a little, 
we can write 
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We have seen the first integral. The second gives . 031075 7695/167r^ and the last gives 
.0142016/1677^; so C2 = .0334416. 

Next we must compute the 0{p^) part of the setting sun diagram. The continuum 
diagram is IR and UV log divergent, while the lattice diagram is only IR log divergent. It is 
convenient to IR regulate both by introducing a mass on one line. In this case the continuum 
integral can be performed in MS, leaving a lattice integral minus an analytically determined 
counterterm [24,25]. Choosing to separate the renormalization dependence along with the 
same finite constant as in the previous literature [24,25], the constant C3 is given by 
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The problem with this expression is the logarithm. To remove it, we add and subtract 
Ic{k) = 1/8A; to hik). The integral 
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is IR convergent, and the m — > limit may be taken immediately. It evaluates to 
— .06858432/1677^, unless we use the unimproved lattice Laplacian, in which case it is 
.60953343/1677^. We re-arrange the remaining terms to be 
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k,BZ V A;2 + m2 k'^ + m? 8k Jk,BZ 8k{k^ + 
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Again, for the first integral the m ^ limit may be taken immediately, and the numerical 
value is .16 1 799 6 07/167r^, or .43364112015/167r^ if we use the unimproved Laplacian. For 
the last integral, we cut the integration region into the ball of radius 7r and the region within 
the Brillouin zone but outside the ball: 
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The first integral is easy and gives ln(7r/m)/167r^ plus terms power suppressed in m. When 
added to (-l/167r2)(ln(6/m) + 1/2), this cancels the ln(m), leaving (l/167r2)(ln(7r/6) - 1/2). 
The final integral has had the small k part of the integration range removed, so the m — > 
limit is trivial. It can then be reduced to 

1 f f^Q 1 12 /"'r/4 /■arctan(sec 0) 

-— / — ln(i?(cube) - i?(ball)) = r-r^— / dcf) / sin(^) ci^ln(sec(^)) (B33) 

which numerically equals .19233513195/167r^. Note that at no point have we had to deal 
numerically with an integral which is log divergent in m, or which still contains m at all. 

Combining terms gives C3 = —.86147916, unless we use the unimproved lattice Laplacian, 
in which case it is C3 u = .08848010. 



4. Three-loop result for S^^ 

In principle, one could extend all the renormalization counterterms we have considered 
to three loop order. However, for many of them, this requires including mixing of different 
dimensions of operator insertions and including counterterms for radiatively induced high- 
dimension operators in the Lagrangian. However, the three-loop contribution to the additive 
countertcrm is an exception: it gives the 0(a) corrections to A(0^)c, and we've already 
seen that we can ignore higher- dimensional operators at this order. The calculation of 
is the least complicated of the 3-loop calculations we might envision, and does not require 
the result of any other 3-loop calculations or any of the new counterterms which would be 
needed in a complete 3-loop matching. The relevant diagrams are given in Fig. 19. We find 
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Here the new constant C4 arises from the basketball diagram of Figure 19, together with part 
of the mass-squared counterterm diagram. The explicit renormalization scale dependence 
in (B34) cancels the implicit dependence of r = r(/2) in the last term. 

Our definition of C4 is that (C4+.^ log(a/i))/ (47r)^ equals the lattice value of the basketball 
diagram, (A) in Fig. 19, minus diagram (C) taking only the setting sun part of the two loop 
mass counterterm. Explicitly, 



C4 + Clog(a/^) 
(47r)3 




p,BZ p2 
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Mp) 



2 HP+k) . 



(B35) 



Here Ii^{p+k) and Ic{p+k) arise from the basketball diagram, while Ih{p) and Ic{p) are from 
the setting sun piece of the counterterm. Every term here is implicitly IR regulated by a 
common infinitesimal mass on every propagator, but the way we will perform the integrations 
means that we will never need this regulation explicitly. Also, continuum integrals are 
implicitly renormalized in MS, which will be relevant. 

It is convenient to add and subtract the appropriate factor to put the Ic{p+k) factor 
inside the BZ k integral. We add and subtract 
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The p integral can be evaluated in MS, (note that we do not need its value in the deep IR 
where the implicit mass regularization becomes important), and gives 
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1 , 1.5 - log(A;//x) 
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p,si3 p^ IGtt^ 



(B37) 



The constant terms can be pulled out of the k integrations; the remaining k integral is 
Eq. (B6) and gives ^/47r. The constant parts therefore yield (1.5 + log(a/i))^/(647r^). The 
In k leads to the integral 




FIG. 19. Three-loop vacuum diagrams needed for Scp"^ at 3 foops. There are 7 additional diagrams either 
involving 1 loop mass counterterms or tadpoles, which cancel among themselves. 
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-ln(aA;) _ r -ln(afe) \ ^ 0.30837 
167r2 \^yfc,BZ A,Si3 J 647r3 ^ 

for the improved Laplacian. Eq. (B35), after subtracting Eq. (B36), becomes 

It is convenient to split off tfie part of the continuum p integration which hes outside the 
Brillouin zone as a separate integration, which does not suffer from divergences: 



-/ 



.00031757 , 

(B40) 



,BZ (A;2)2 yp_sR3_Bz p'^\p+k\j 647r^ 

This leaves as the final integral we must consider, 

L (k L (f ^'^^^'^ - '^'"M - ? ^'^^^'^ - '^^") ■ 

This integration is finite and over a finite integration region (provided we make a prescription 
that the argument of the p integral is interpreted as the average over p and —p of the 
quantity written). The integration is nine dimensional and contains some delicate integrable 
singularities. We find it convenient to use adaptive Monte Carlo integration rather than 
quadratures. Monte Carlo integration is somewhat problematical when there are integrable 
singularities. Rather than finding clever changes of variables to get rid of such singularities, 
the simplest thing to do is to include a small mass in all the propagators (thus cutting off 
all singularities), vary that mass, and then numerically extrapolate the zero mass limit from 
the results of the Monte Carlo integrations. Our result for the integral (B41) is (0.0985 ± 
0.006)/(47r)^. Combining terms, we then have C4 = 0.2817(6) for the improved Laplacian. 

5. Minimalist expression for A(0^)c through 0(a) 

Some of the development in the preceding sections is unnecessary for the isolated goal 
of getting an 0{a) improvement of A(0^). It is possible to combine the previous results in 
the more compact form 

C = C4 + ^ln6 + ^C3. (B42b) 

Here $ is the bare lattice field, and Uq the bare lattice coupling, corresponding to the bare 
lattice action (2.1); 

NT. 

A($^).at = ($')lat - 1^ ; (B42C) 
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ri stands for the tadpole-adjusted mass, 

n^r,+ ^^^u,^^ (B42d) 
Ana 

which is 0{u^) near the transition; and tq is the bare mass used in the bare lattice action 
(2.1). We have retreated from lattice units and explicitly show all factors of a. Since A($^) 
is 0{u) near the transition, we only need the one- loop result 

(B42e) 

for Zr/Z^. 

To finish the relationship between the measurement of A(0^)/m and the bare fields and 
parameters of the lattice Lagrangian, we only need the one-loop relationship between the 
continuum and bare lattice couplings: 

u^^u+^^^^u^a^. (B42f) 
6 An 

The form (B42) has the conceptual advantage of not introducing the renormalization 
scale Jjb, since its introduction is unnecessary if one's only interest is in A(0^)c and not the 
value of Tc- Eq. (B42) also makes clear that nothing depends on Ci and C2 at this order, 
so that one doesn't really need those integrals. And it makes clearer that there is no 0(a°) 
correction to A(0^), which is obscured by eqs. (2.15) and (2.16). That is, 

A(02)eont = A($2)iat + 0{au''), (B43) 

as in (2.11). 

However, this is not how we have implemented our calculation of A(0^)c when quoting 
numbers from simulations. What we have done is described in the text and the earlier parts 
of this appendix and, though equivalent through 0(a), will give slightly different numerical 
values because of differences in higher orders in a. 



APPENDIX C: LOGS IN LARGE VOLUME SCALING FOR a = 

In this section, we review standard renormalization group arguments about the free 
energy and verify that a In L term appears in the free energy when a = 0. If one increases 
renormalization scale by a "blocking" factor of 6, the free energy density of the blocked 
system is related to the free energy of the original system by a transformation law of the 
form 

r{{K}) = g{{K}) + h-''H{K']l (ci) 

where T is the free energy per block, {K} represents the couplings of the theory, and Q is 
an analytic function (depending on h) that represents the contributions to the free energy 
from the degrees of freedom that have been blocked.^ Iterating n times, this becomes 



^Our presentation follows, for example, chapter 3 of Ref. [26], which is one of many nice introductions to 
the renormalization group. Our T and G are that reference's / and g. We use T to avoid confusion with / 
in the text, which was the free energy per physical volume rather than per block. 
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n-1 

= E b-"'gi{K^^'^}) + (C2) 

j=0 

where {K^^^} is the jth iterate of {K}. If we start with a block size of a and iterate all the 
way out to size L, we have n — logj,(L/a). 

For t very small but non-zero, the description of the system will first flow towards the 
critical point, as one blocks to larger and larger distances. But it will then eventually flow 
away from the fixed point, closely following one of the two unique "outfiow" trajectories from 
the fixed point (one for t > and one for t < 0), corresponding to all irrelevant couplings 
being set to zero. The inhomogeneous part Q of the transformation law (C2) will generate 
singular behavior as i — > 0, which may then be written as 

logj (L/a) 

^(t) ~ E b-^'^Giy^H), (C3) 

j=0 

where IP'^^t parameterizes fiow along the outfiow trajectory. For the infrared behavior at large 
L and small t, we can replace the sum by an integration. Changing integration variable to 

s = vy% 

fd/yt i-t(L/a)yt 

m - ^ / s^~'g{s) ds, (C4) 

where a — 2 — d/yt- For a — 0, this is 

To help tame the singularity as s — > 0, integrate by parts two times. This leaves an integral 
proportional to 



s-^g"{s) ds (C6) 



plus terms which fit the too-naive scaling form (4.12). The remaining integral can be rewrit- 
ten as 

Q"{^) ln(L-'^/2) + [g"{s) - g"{0)] ds. (C7) 

The second term has the desired analyticity properties of the too-naive scahng form (4.12), 
and the first term is the logarithm of Privman and Rudnick appearing in (4.14). 



APPENDIX D: OUR OWN ANALYSIS OF Cc 

1. Numerical simulation 

In this appendix, we discuss our own attempt to determine Cc, to make a crude check 
of the value we take from Ref. [13]. The scaling of the intersection values C^{Li,L2) of 
cumulant curves C{r) is 
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fit 




C.L. 


Ls/a > 16 


1.229(9) 


NA 


Ls/a > 12 


1.243(4) 


6% 


Ls/a > 8 


1.2404(22) 


13% 


Ls/a > 6 


1.2410(15) 


25% 


Ls/a > 4 


1.2402(7) 


35% 


Ls/a > 3 


1.2394(7) 


1% 



TABLE I. Results of linear fits to the subset of data of fig. 20 with Li = 2L2, and the corresponding 
confidence levels ("NA" means not applicable). The first row corresponds to a fit to two points, the second 
to three points, etc. 



j-yt r -uj j-yt T-uj 1-u) vyt 



for Li, L2 00. This is a simple consequence of Binder's analysis [9], but, for the sake of 
completeness, we briefly outline the argument in Appendix D 2. 

We would Ukc to choose L as large as possible, in order to make the scaling law (Dl) as 
accurate as possible, so that we can use it to extrapolate a good value of Cc- Our ultimate 
interest in this paper is to study continuum 0(2) theory, which requires ua <^ 1 and for 
which the scaling limit is L/a ^ 1/ua. But L/a':$> 1/ua implies that small ua simulations 
are an inefficient choice for getting as far as possible into the scaling limit. Because Cc is 
universal, we can extract it from large ua simulations rather than small ua simulations. [In 
fact we could use any model in the same universality class.] Our simulation code is optimized 
to perform best if ua is not extremely large, and so we have chosen to extract Cc from data 
taken with uoa — 60. Because uoa — 60 is not a small value of ua, the 0{a) improvements to 
the action are pointless; unlike other simulations reported in this paper, we quote quantities 
in terms of the bare lattice parameters. 

A variety of intersection values Cx{Li,L2) of C(r) curves are plotted in Fig. 20 against 
{Ll'L2^ — L2*i>f'^)/(Lf — L2*) which, by (Dl), should lead to a hnear relationship at large 
LxiLi- The errors on points sharing an L value are correlated, and we compute and use the 
full correlation matrix for making fits. For the sake of simplicity, however, wc have only fit 
the subset of data with L\ = 2L2, with results for Cc and the associated confidence levels 
given in Table 1. If we naively include smaller and smaller L until the confidence level of 
the fit becomes poor, we would fit all the way down to L = 4 and obtain 

Cc = 1.2402(7) (D2) 

as our final result. We are suspicious that this result may have a "tail wagging the dog" 
systematic error. Our small L data has smaller statistical error than our large L data and 
may be over- weighted in the fit to what's supposed to be large L asymptotic behavior. (In 
particular, in an earlier analysis where we had even poorer statistics on the large L data, 
we found a result that disagreed even more drastically from the result of Ref. [13].) 
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FIG. 20. The values Cx{Li, L2) corresponding to the intersection points of pairs of the curves C(r) 
for different lattice sizes. The horizontal axis has been chosen so that this relationship should be linear as 
Li, L2 00 (which corresponds to approaching zero on the horizontal axis). Each point is labeled by Li/a 
and ^2/1. The linear fit shown is only to that subset of points with Li = 2L2. 



2. Large volume scaling of cumulant intersections 

Let's quickly reproduce Binder's result (3.5) for the scaling of rx(Li, L2) — r^. For large 
L, we will focus on the scaling piece of the cumulant which, in the notation of section IV A, 
is 



Choosing b = L, 



C{t, {u,}, L-i) ~ C,i„g(6^'t, b/L). 



C~Csing(L^%KL%},l), 



(D3) 



(D4) 



where Cging is a universal scaling function. Now treat L^^t and UjL^^ as small and Taylor 
expand, keeping track of only the most important irrelevant operator uj, 



C ~ aing(0, 0, 1) + AUH + BL-'^. 



(D5) 



Note that A and Csing(0,0, 1) are universal, but B is not, because it depends on the value 
of Uj. Now look for the intersection for two different system sizes: 



aing(0, 0, 1) + AL\H + BL^^ = aing(0, 0, 1) + ALlH + BL:,^, 
which has solution 

B{Lr - 



t 



A{Lf -Lf 



(D6) 



(D7) 



up to corrections suppressed by additional powers of 1/L. This is just the scaling (3.5) 
quoted for rx — Vc- (And one sees a posteriori that it was justified to treat L^H as small in 
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the L ^ oo limit.) All that is necessary to derive the scaling (Dl) of Cx is to plug (D7) 
back into (D6) and rename Csing(0,0, 1) as Cc- 

In the main text, our procedure for defining the nominal transition in finite volume, for 
our small ua simulations, is to find the point where C = Cc- Prom (D6), one may verify 
that this prescription gives 

t ~ ^ L"^*-^ ~ L"^*-^, (D8) 

as asserted in the main text. 



APPENDIX E: SMALL Lu EXPANSION OF A(02)^ 

In this appendix, we compute the expected result for A(</)^)c for small volumes in the 
continuum limit, where we define the nominal critical temperature in such volumes using 
Binder cumulants, as explained in Sec. III. We have repeatedly described "short distance" 
physics (distance small compared to 1/u) as perturbative, and so for small volumes (L <S 
1/u), one might first try simple perturbation theory to compute A(0^). The leading-order 
diagram contributing to A(0^) is shown in Fig. 16d and corresponds to 

where for the sake of generality we have considered an O(A^) model, and the actual case of 
interest is N — 2. The sum is over the discrete momenta p — 27rn/L associated with the 
Lx Lx L periodic volume. There is a problem with (El), however. Consider the mean- field 
theory approximation Tc = to the critical value of r. The sum then has a infrared divergent 
term associated with p = 0. The problem is that it is more accurately large momentum 
physics which is perturbative, rather than small distance physics. So, even when L — > 0, the 
physics associated with the p = modes is non-perturbative. Since this is only one mode, 
we simply analyze it separately, and then treat all the p 7^ modes perturbatively. 



1. The case r — 0: leading order 



To illustrate the calculation, let us first compute A(0^) if r = 0. We'll later come back 
to consider the actual r that is chosen by the criteria of our simulations that the cumulant 
C — Cc- Consider the approximation So to the action where we ignore everything but the 
p=0 modes 



drx 



(E2) 



Note that this result holds on the lattice as well as in the continuum. The p = contribution 
to A(02) is then 



So 



Lht 



1/2 



/ d X e ^ X 



Lht 



r( 



JV+2 



r(f) ■ 



(E3) 
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By dimensional analysis, the leading perturbative contribution of the p 7^ modes to {(fP') 
should be order L^^ and so is dominated by the zero-mode contribution (E3). Specializing 
to iV = 2, we then have 

As we'll see explicitly below, the UV subtraction that converts (0^) to A(0^) in the contin- 
uum limit is not relevant until next order in (Lm)^^, and so A(0^)/m also has the expansion 
(E4). 



2. The case r = 0: next-to-leading order 

Now consider the leading-order perturbative contribution 5{<p'^) of the non-zero modes 



to (02) from (El): 



5^ 



V — 



(E5) 



This sum is UV divergent, and the UV divergence is simply the free-field divergence discussed 
in section (HA). It is subtracted when we compute A(02) as opposed to (0^). Note that 
our prescriptions (2.11) or (2.15) for computing A(02) in our simulations involve subtracting 
the infinite volume free-field result for (0^). In the continuum limit, this subtraction turns 
the perturbation (E5) into 



5 A 



N 
1? 



(E6) 



Individually, the sums and integrals above must be consistently regulated, for example by 
dimensional regularization or by keeping the system on a lattice with arbitrarily small lattice 
spacing. 

There are a number of ways to evaluate the result numerically. One is to start by 
regulating with dimensional regularization, working in d spatial dimensions. Then 



''A(0^> = ^Ei-A'/i 



N 



/pp2 (27r)2L'^ 



-2 



?n^' / 



(E7) 



Now rewrite 



Jo 



ds e 



The ni, ^2, ... sums and integrals then factor, giving 



Y.-- 



d'^n 



ds 



- 1 - 



7ry/2~ 



(E8) 



(E9) 
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where 

oo 
fc=— oo 

is a special case of an elliptic theta function. The integral on the right-hand side of (E9) is 
absolutely convergent in d = 3, so we can now dispense with dimensional regularization and 
set d to three. The integral is then easily done numerically, giving 

j^l_fd^^ -8.91363. (Ell) 



For N — 2, one then has 



S^{<t>')r=o ^ 0-451570 
u Lu 

The effect of interactions in the non-zero mode sector will be perturbative for small 
Lu and will give higher order contributions to A(0^). However, we must also consider 
interactions of the zero mode 0o with the non-zero modes, which will give corrections to 
the action 5*0 used in our earlier analysis of 0o- Foi' instance, the coupling u in 5*0 will pick 
up corrections of order Lv?. This will generate an 0{Lu) relative correction to the zero- 
mode contribution (E4), and so that correction will be higher order than the contribution 
(E12) computed above. A similar story will hold for corrections to r (and for the effects 
of higher- dimensional interactions) with the complication that r will receive some infinite 
contributions in the continuum limit, corresponding to the usual mass renormalization. The 
latter is absorbed by the usual renormalization of r, and so one has 

.(^^ 276m _ 0451570 ^ 

u [Luyi'^ Lu 

if one interprets the r in the condition r = as being renormalized r. (The details of the 
renormalization scheme will not matter at the order shown.) 



3. The case C — Cc 

Now, instead of setting r = 0, we will choose r so that the cumulant C is equal to its 
critical value. Let's begin with a leading-order analysis and so focus on just the p = sector. 
For general r, we then have 

S.-^\<t>o? + ^\M' (E14) 

and 

_ /ci^0oe-^o|0oP^ 
By a change of variables, this can be rewritten as 
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where 



and 



(A\\"'h{R) 



UR) = r dy e-2/'-«^y'=+(^-2)/2^ (E17) 

^0 



) . (E18) 

Now note that the volume average of the field, used in the definition of the cumulant, 
is simply 0o- Specializing now to A?" = 2 and doing the 7^ integrals, one obtains 



2 V2 

h{R)^\[l-RIo{R)l (E20) 
h{R) ^\[-R+ {R^ + 2)Io{R)]. (E21) 



The cumulant can be written 



(4) ^, h{R)m) 



A numerical search for the point where C — Cc — 1.243 gives 

i?~ -2.5073. (E23) 
At this R, the zero-mode contribution to (0^) is then 

Since the result for i? is a pure number, the definition (E18) of R shows that r is of 



order ^u/L^, which is smaller by a power of y Lu than any non-zero momentum squared, 
which are order L~^. So r can be ignored in finding the leading contribution of the non-zero 
modes, with the effect that 5A(0^) is the same as in the earher r = analysis. Adding 
(E24) and (E12) then produces the result (4.15) quoted in the main text. 



APPENDIX F: SYSTEM SIZE AND THE SIMULATIONS OF GRUTER et al. 

One of the applications of 0(2) field theory is to studying the corrections to the criti- 
cal temperature, due to interactions, for Bose-Einstein condensation of a nearly-ideal non- 
relativistic Bose gas. This application has been previously studied using numerical tech- 
niques. Our discussion of system size in section IV C has implications for an early study 
by Griiter et al. [22], which is that much of the data collected was likely in insufficiently 
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large volume. These simulations did not make use of 0(2) field theory: They worked in the 
canonical ensemble and studied the path integral for a fixed number A^p of particles in a finite 
volume V. If no attempt were made to fit for large volume corrections, then Fig. 2 makes 
it clear that one should take roughly Lu > 400 to keep those corrections moderately small. 
It's illuminating to also consider the less restrictive condition of roughly Lu > 200. We can 
translate these conditions on system size to the context of Bose-Einstein condensation by 
translating the parameter u of the 0(2) field theory. For this application, the relation is 
[2] u = 967r^a/A^, where a is the scattering length of the atoms, A = hJ2n/mkBT is the 



thermal wavelength, and m is the mass of the atoms. Using the ideal gas approximation 
Tc ~ To = (27r/i^//cBr«)[n/C(|)]^/^ for the critical temperature, one can write 



at the transition. The conditions Lu > 400 or 200 may then be translated into the conditions 

o 0.51 0.064 

The largest Np used in the simulations of Griiter et al. was A^p = 216, and their extraction of 
the critical temperature depended on results with A'p = 125 as well. For the latter, our rough 
conditions for being in large enough volume then requires na^ > 0.004 or 0.0005, depending 
on whether one takes the more or less restrictive condition Lu > 400 or Lu > 200. Griiter et 
al. quote results for na^ all the way down to 10~^, and the majority of points in their small 
na^ tail have na^ < 10"*^. (See Fig. 3 of Ref. [22].) It therefore seems at least possible that 
they may have had inadequate system sizes for their extrapolation of the — > behavior. 
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